%% 
clear;
close all;

%% 定义信道参数
%--------------------------------------------------------------------------

M_mod = 4;                         % 4QAM
k = log2(M_mod);                    % Bits/symbol
M = 512;
Mzp = 180;
B = 5000;                      % 带宽 Hz
scs = 10;                      % 子载波间隔 Hz
ofdmSym = 16;                   % ofdm信号数量 /帧 
fc = 12500;
c = 1500;                                                     % speed of light (m/s)


%% channel 2
N = ofdmSym;

Fn = dftmtx(N);       % Generate the DFT matrix
Fn = Fn./norm(Fn);    % normalize the DFT matrix

Fn2 = dftmtx(M);       % Generate the DFT matrix
Fn2 = Fn2./norm(Fn2);    % normalize the DFT matrix

Md = M-Mzp;

data_grid = zeros(M,N);
data_grid(1:Md,:) = 1;
N_syms_perfram = sum(sum(data_grid));

delta_f = scs;
T = 1/delta_f;
% length_CP=cpLen;

one_delay_tap = 1/(M*delta_f);
one_doppler_tap = 1/(N*T);


max_speed=20;  % km/hr
[chan_coef,delay_taps,Doppler_taps,taps]=Gen_delay_Doppler_channel_parameters(N,M,fc,delta_f,T,max_speed,7);
L_set=unique(delay_taps);
[G,gs]=Gen_time_domain_channel(N,M,taps,delay_taps,Doppler_taps,chan_coef);



             
%%                       OTFS BER Calculation

     data_na= randi([0,1],N_syms_perfram*k,1);
     data=qammod(reshape(data_na,k,Md*N), M_mod,'gray','InputType','bit');
     X = Gen_2D_data_grid(N,M,data,data_grid);
    
    
    %% OTFS modulation%%%%
        X1_tilda=X*Fn';
        s1= reshape(X1_tilda,N*M,1);
    
     %% channel output%%%%%

% %     r1=zeros(N*M,1);
%     noise1= sqrt(0.2/2)*(randn(size(s1)) + 1i*randn(size(s1)));
%     l_max=max(delay_taps);
%     for q=1:N*M
%         for l=(L_set+1)
%             if(q>=l)
%                 r1(q)=r1(q)+gs(l,q)*s1(q-l+1);
%             end
%         end
%     end
    r1=G*s1;
    [r2,tau_p,A_p,a_vel] = MChannelchirp2(s1,5,fc,B);
    r2=r2(1,1:M*N)';

    delay_taps2 = round(tau_p/one_delay_tap);
    Doppler_taps2 = a_vel/one_doppler_tap;
    L_set2=unique(delay_taps2);
    [G2,gs2]=Gen_time_domain_channel(N,M,taps,delay_taps2,Doppler_taps2,A_p);
    r22=G2*s1;

%     figure(1)
%        plot(0:1/B:(M*N-1)/B,abs(r2),'y',0:1/B:(M*N-1)/B,abs(r22),'r');
%     legend('MChannelchirp方法生成信道','使用MChannelchirp参数-OTFS方法生成信道')
%     title('输出信号')

%      figure(2)
%      [R,log]=xcorr(abs(r1),abs(r22));
%      plot(log,R/M*N);
%      title('输出互相关')

    

    figure(3)
    plot(0:1/B:(M*N-1)/B,abs(r1),'y',0:1/B:(M*N-1)/B,abs(r2),'r');
    legend('MChannelchirp方法生成信道','OTFS方法生成信道')
    xlim([0,0.4]);
    title('输出信号')
% 
    figure(4)
    [R,log]=xcorr(abs(r1),abs(r2));
    plot(log,R/M*N);
    title('输出互相关')
 
  
     %% OTFS demodulation%%%%
    Y1_tilda=reshape(r1,M,N);
    Y1 = Y1_tilda*Fn;
    y1=reshape(Y1.',N*M,1);

    %% detector compare%%%%
    eng_sqrt = (M_mod==2)+(M_mod~=2)*sqrt((M_mod-1)/6*(2^2));
    SNR_dB = 5:5:25;
    SNR = 10.^(SNR_dB/10);
    noise_var_sqrt = sqrt(1./SNR);
    sigma_2 = abs(sqrt(eng_sqrt)*noise_var_sqrt).^2;     

    errorRate1 = comm.ErrorRate('ResetInputPort',true);
    errorRate2 = comm.ErrorRate('ResetInputPort',true);

    berOTFS = zeros(length(SNR),3); 
    berMChannel = zeros(length(SNR),3);

    for m = 1:length(SNR)
        noise= sqrt(sigma_2(m)/2)*(randn(size(s1)) + 1i*randn(size(s1)));
        r1_n=r1+noise;
        r2_n=r2+noise;
        r22_n=r22+noise;

        [est_LMMSE_r1_n,~] = LMMSE_detector2(N,M,M_mod,sigma_2(m),data_grid,r1_n,G); 
%         [est_LMMSE_r2_n,~] = LMMSE_detector(N,M,M_mod,sigma_2(m),data_grid,r2_n,gs2,L_set2); 
        
        berOTFS(m,:)=errorRate1(data_na,est_LMMSE_r1_n,0).';
%         berMChannel(m,:)=errorRate2(data_na,est_LMMSE_r2_n,0).';
    end
%% figure
berOTFS(~berOTFS)=1e-12;
% berMChannel(~berMChannel)=1e-12;

figure
hold on;
semilogy(SNR,berOTFS(:,1),'-r');             %Plot simulated BER w/ OTFS
% semilogy(SNR,berMChannel(:,1),'--b');             %Plot simulated BER w/ C-OTFS
    

ylabel('BER');
xlabel('SNR');
legend( 'berOTFS');
title("输出误码率")
grid on;
hold off;


